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ABSTRACT 



Subspace methods of solving spectral estimation and direction of arrival (DOA) 
problems involve finding the eigenvalues and eigenvectors of correlation matrices. Using 
the Lanczos algorithm some of the extreme eigenvalues and eigenvectors can be ap- 
proximated without requiring the entire matrix decomposition theoretically saving many 
computations. 

This thesis develops a model and a form of the Lanczos algorithm to solve the DOA 
problem. The relationship of the number of eigenvectors used to the accuracy of the 
results in a low signal to noise ratio example are examined. 
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I. INTRODUCTION 



Recently, a number of signal processing techniques have been developed that in- 
volve resolving a received signal into so-called signal and noise subspaces. The ability 
to perform spectral estimation or direction of arrival from the determination of the noise 
subspace has been shown in the works of Pisarenko, Schmidt, Kay, Kailath, and others. 
The methods vary in the manner in which the subspace is reached and how the resulting 
eigenvectors are calculated. [References 1,2, and 3| 

To achieve better results (detection at lower signal-to-noise ratios, better resolution, 
finer accuracy, less bias), more samples are required. This leads to larger, more complex 
matrices that better describe the signal and noise subspaces. Determination of the sub- 
spaces requires an eigendecomposition of an estimated correlation matrix. The general 
eigendecomposition of a matrix requires computations 0{r&\ thus the larger matrices 
require many more operations. Once the matrix is decomposed into its eigenvalues and 
eigenvectors the proper set of eigenvectors must be used to find the resulting spectrum. 
Hence, estimation is required to split the eigenvalues into signal and noise subsets. 

The difficulties in the procedures can be stated with two questions. 

• Where is the threshold between signal and noise eigenvalues (and thus which, or 
how many eigenvectors are used)? 

• What method should be used to find those eigenvectors in a timely fashion? 



Proposed here is a method which will accurately estimate some of the eigenvalues 
and eigenvectors of a matrix without performing the entire decomposition. Computa- 
tional savings are realized when only a partial decomposition is required. The 
eigenvectors used correspond to the minimum eigenvalues of the matrix. With an 
over-specified matrix (dimension much larger than the number of signals present), these 
minimum eigenvalues will all correspond to the noise subspace. Each noise eigenvector 
contains all the information to find the actual spectrum, although spurious results will 
also be apparent (because the matrix is over specified). Several methods of handling 
these spurious peaks are illustrated, including eigenvector averaging, spectral averaging 
and using the spectral product. 

This thesis is organized in five chapters. Chapter II, Direction of Arrival, discusses 
classical spectral estimation theory and how it applies to the linear array problem 
(beamforming). Subspace methods starting with Pisarenko Harmonic Decomposition 
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and proceeding to MUSIC and ESPRIT are discussed in detail. Chapter III, The 
Lanczos Algorithm, includes all the basic theory required to describe this 
eigendecomposition method. There we also compare several methods to negate the 
spurious peaks with the proposed spectral product giving the best results. Results of the 
algorithm for numerous cases are given in Chapter IV. The last chapter summarizes the 
results and advantages of this Lanczos subspace method. This chapter also includes 
some recommendations for future work and lists some possible applications. 
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II. DIRECTION OF ARRIVAL 



Direction of arrival is a form of spectral analysis, performing frequency detection 
and resolution in the spatial domain vice the conventionally considered time domain. 
The signals incident on an array are analyzed, and, if the presumptions of the analysis 
are valid, the correct bearings to the sources are determined. Formerly it was only pos- 
sible to analyze the output of an array by conventional Fourier techniques. More re- 
cently numerous methods have been developed which enhance the ability to accurately 
determine the spectral and angular resolution. 

This chapter summarizes some of the salient points of spectral estimation before 
discussing classical direction of arrival array processing (Bartlett beamforming). 
Projection-type superresolution subspace methods are then discussed, starting with 
Pisarenkos Harmonic Decomposition. MUSIC and ESPRIT are discussed in detail and 
several other subspace methods are mentioned. 

A. SPECTRAL ESTIMATION 

Spectral estimation is the term used to describe efforts to find the frequency com- 
ponents of a signal sampled in time. Two conditions that are required for the remainder 
of this thesis are that the processes considered are wide sense stationary (WSS) and 
ergodic. The assumption of WSS means that the process has finite mean and that the 
autocorrelation is a function of the distance, or lag, between two samples and not of the 
position of the samples themselves. Ergodicity allows time averages to be used to de- 
termine ensemble averages. 

The autocorrelation function of the signal x(t) is 



where jc(m) are the individual samples of the signal. When only a finite number of sam- 
ples, M, are taken, the above definition must be modified. The sample autocorrelation 
function is defined as 




(1) 



3 



( 2 ) 



*«(*) 



M- 1- 

■+z 



jc(« + ife)x(rt), for 0 < & < (A/ — 1) 



It is easily shown that [Ref. 4: pp. 56-58] 



*«(*) = * kc (-*). for -( M - 1)<*<0 
and 

4,(0) > 4*w. for aii * 



( 3 ) 



The sample autocorrelation matrix R„ is 



R xx (0) R xx {- 1) 

4,d) 4,(0) 



RJW- 1) 



R„(-A/+l) 



*„( 0 ) 



(4) 



The power spectral density (PSD) is a measure of power per unit frequency. A plot 
of PSD versus frequency will show the relative power at all the frequencies present. It 
also describes the properties of the noise in the signal, i.e., white noise will have a flat 
PSD showing that all frequencies are equally represented [Ref. 1: pp. 53]. The PSD is 
given by 



(/) - Y Yj R ^ e ~ J2n/ " T 



where T is the sampling period. 

The periodogram method of estimating the true PSD is one of the earliest and most 
widely used [Ref. 5: pp. 5-8]. The periodogram is defined as the z-transform of the 
sample autocorrelation matrix evaluated on the unit circle [Ref. 4; pp.52-53] 
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M - 1 

s„(z) - ^ R„(k)z- 

k=-(M- 1) 



(6) 



It may also be found by directly z-transforming the original data sequence x(n) 



M 

$xx( z ) - -yf X{ z )X(z~ l ) where X(z) = ^ x{n )z~ n 



( 7 ) 



The periodogram spectrum is found by substituting z = e 12nfT , 



L(f) - -^U'WI 2 = -L 



M 

X 



x{n)e J 



(S) 



Data is often run through a computationally efficient Fast Fourier Transform (FFT) to 
find the periodogram spectrum. 

The measures of effectiveness of a spectral estimator are based on comparisons of 
resolution, detectability, bias and variance. Resolution relates to the ability of the esti- 
mator to separate two separate spectral lines that are closely spaced in frequency. The 
capacity to locate a low energy signal is measured in detectibility. Bias is the tendency 
of the observed peak to be shifted off the true spectral line. Variance is a measure of the 
propensity to keep the true spectral shape, or how quickly the PSD falls off away from 
the true peak. Different spectral estimators are sometimes compared in terms of 
robustness, or ability to function well with various types of signals and noise. 

If individual signals are narrowband the resolution criterion is said to be achieved 
when direct observation of the spectrum leads to the correct determination of the num- 
ber of signals. Separate peaks are not required to determine that two signals are present 
[Ref 6]. Resolution is inversely proportional to the length of the data samples. With f 
as the sampling rate and T the sample period, or time between samples, T = -jr, the fre- 
quency resolution is given by 
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A r _ Jj_ _ _!_ 
^ M MT 



(9) 



Thus, with the periodogram, better resolution can only be achieved with longer record 
lengths. 

The above description is based on no windowing (or rectangular windowing) of the 
data samples. The windowing in the time domain is seen as convolution in the frequency 
domain. The rectangular windows, for example, transform into sine functions in the 
frequency domain resulting in high sidelobe distortion known as leakage in the frequency 
domain. High sidelobes result in many false peaks, making it difficult to discern the true 
spectral peaks, so detectibility suffers. 

Nonrectangular windows are used to taper the data to minimize the discontinuities 
that cause the high sidelobes. Common windows include the Bartlett and Hamming. 
The lower sidelobes come at a cost of resolution, so two or more signals close to each 
other in frequency may merge into one. Resolution may be boosted by lengthening the 
sample time, but the increased record length may violate the requirement for 
stationarity. More can be found on the subject of windows in References 7 and 8. 

Because of the difficulty in meeting the requirements for both detectibility and re- 
solution, parametric methods have been developed that produce increased resolution 
with short data lengths. The parametric methods are based on determining an appro- 
priate model for the process that produced the data. If the process can be effectively 
modeled, then more reasonable assumptions may be made about the data's behavior 
outside of the window (i.e., these data points do not have to be set to zero). Once the 
appropriate model is chosen to represent the process, the parameters are estimated from 
the available data, inserted into the model, and the power spectral density expression for 
the respective model is determined. A few common parametric methods are summarized 
below. [Ref 1: pp. 106-108, Ref 5: pp. 172-174] 

Many random processes that are encountered in practice may be represented by the 
linear difference equation 



p 



<7 




(10) 
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where u{ri) is the input driving sequence and x(n) is the output sequence. The a(k ) pa- 
rameters are the autoregressive coefficients and the b(k) are the moving average coeffi- 
cients. Equation 10 is thus known as the autoregressive-moving average (ARMA) model 
or ARMA r p,q) process and is the most general model with a rational transfer function. 
The power spectral density that results from the ARMA model exhibits both sharp 
peaks and deep nulls. If the autoregressive parameters of equation 10 are all set to zero 
except a(0) = 1, the resultant model is the moving average (MA) process of order q. The 
MA spectrum will have deep nulls, but relatively broad peaks. With the b(k) coeffi- 
cients of equation 10 set to zero except b(0) = /, the autoregressive (AR) process of order 
p results. The AR spectrum will exhibit sharp peaks but will have broad nulls. [Ref 5: 
pp. 174-181] Each of the models will exhibit greater accuracy and flexibility as the order 
is increased. With a high enough order the AR method can approximate an ARMA or 
MA process, and, likewise, a very high order MA model can be used to approximate an 
ARMA or AR process. But if the model order is too high for the actual process, esti- 
mation errors will lead to nonzero coefficients and spurious peaks will result. Thus 
proper estimation of model order is important [Ref 1: pp. 112-113, pp. 19S-207]. 

The parameters of these three models may be estimated by using the Yule- Walker 
equations utilizing the autocorrelation matrix of the available data stream [Ref 1: pp. 
1 15-1 IS] 



R^a = -r (11) 

While the true autocorrelation function is impossible to determine, the maximum 
likelihood estimator (ML) is one means to find good approximations of the parameters 
for the AR model. The ML method uses a suitable estimate of the autocorrelation or 
covariance matrix and then solves [Ref 1: pp. 185-190] 

Ca — — c (12) 

for the parameters a where C is the symmetric covariance matrix and c is an 
autocorrelated vector. 

The Burg method (maximum entropy) estimates reflection coefficients from the 
process and then uses the Levinson recursion to find the estimated parameters [Ref 1: 
pp. 228-231]. 

Generally, the various AR spectral estimators except the autocorrelation method are 
unbiased and have similar variance. The covariance and Burg methods are restricted to 
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ranges that keep the summation within the available data and do not assume zero pads 
outside of the samples, thus taking advantage of the basis which led to the creation of 
the parametric methods in the first place. [Ref 1: pp. 240-252] 

B. BEAMFORMING 

A conventional approach to the direction of arrival (DOA) problem is through the 
classical beamforming (Bartlett) technique. Basically, this is a measure of coherency of 
the signals arriving at an array of sensors. The characteristics of an array are described 
in terms of array gain, directivity, resolution, beamwidth, and sidelobes. These are based 
on array size (number of sensors), sensor sensitivity, intersensor spacing, and post re- 
ception processing. 

Assuming a narrowband point source in the far field, a plane wave will pass through 
a linear array in a specified order, where the magnitude of the excitation on any indi- 
vidual sensor will be related to the phase of the signal at the instant of sampling. The 
individual sensors will see difTerent instantaneous magnitudes based on this phase dif- 
ference which is a function of the frequency (or wavelength), the DOA, and the spacing 
between sensors. The difference in phase for two successive sensors for a linear array 
with equally spaced sensors is 



where d is the distance between sensors, / is the signal wavelength, and 6 is the angle 
between the wavefront and the array axis. This phase difference A$ is also known as the 
normalized wavenumber, k [Ref. 4: pp. 341-343]. Figure 1 illustrates the array- 
wavefront interaction. 

The output of a simple narrowband delay-and-sum beamformer is 



where jc m (r), m = 1,2,..., M is the signal at the mth sensor. The time lag, r m , between 
two adjacent sensors is to make up for the propagation delay caused by the wavefront 
having to travel the extra distance ^ sin One can adjust the magnitude of the output 
to plane waves arriving from a particular direction 6 by introducing appropriate time 



A , Lna ■ a 
A(p = — : — Sin d 
/. 



( 13 ) 



M 




(14) 



delays at each sensor prior to summing. This is known as "steering the array" or 
"steering the beam". An illustration of a typical beamformcr arrangement is shown in 
Figure 2. 




In the multiple source case, especially if the undesired signals are interfering with the 
detection of other sources, this technique may be modified to steer nulls instead of 
beams thus minimizing output from "jammers". Nulls may be directed toward a single 
known source to minimize the strength of the signal with respect to that source, with the 
output being analyzed to determine what other sources are present. Another imple- 
mentation is to steer the nulls to minimize the total output. The analysis of the delays 
will indicate the directions of multiple targets. [Ref. 4: pp. 351-355, Ref. 9 ] 

Beamforming may also be analyzed in the frequency domain by transforming 
equation 14 into the frequency domain 

M 

«W = xMe-™- (15) 
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In vector notation we can write e. = \v r x where w are the weights and x are the outputs 
of M sensors 









*i(/) 








*2 (/) 


where w = 


. 


and x = 


. 




e ~J2*fr M 




*\M) 



The steering vector s k is the phase relationship between the angle 0 and the normalized 
wavenumber k given in equation 13 
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1 

e* k 



and k = — — sin 6 



(17) 



It can be shown that the steering vector from an array with weights w, directed toward 
an arbitrary’ direction 6, can be expressed in terms of the steering vector as a = s* [Ref. 
4: pp. 343-344], 

Frequency domain beaniforming is directly analogous to spectral estimation 
decribed above. Spatial sampling has the requirement that sensor distances d must be 
less than // 2 apart to prevent "grating lobes" (or aliasing) corresponding to the Nyquist 
rate in the frequency domain of <//2 . Longer arrays, containing more sensors, 
will give better resolution and smoothing. This is similar to frequency resolution being 
proportional to the data record length (A/% ). [Ref. 4: pp. 341-345, Ref. 10: pp. 

4-S, Ref. 1 1 : pp. 293-299] 

The DOA is actually a relative comparison of observed spatial frequency and known 
signal frequency. The spatial frequency is greatest on endfire, when the wavefront is 
perpendicular to the array (0 = 90° or —), Here the phase difference between adjacent 
sensors is at a maximum. A simultaneous sampling of all sensors at one instant in time, 
or snapshot, will indicate the maximum spatial frequency. Observation of the spatial 
wave over time (with a known temporal sample rate) will indicate the end of the array 
where the source is located. 

On broadside (6 = 0 or n), the wavefront excites each sensor identically. Spatial 
sampling indicates no phase difference along the entire array, except for the effects of 
additive noise, resulting in the computation of zero spatial frequency. Unfortunately, 
linear arrays have an inherent ambiguity with broadside signals. The side of the array 
at which the source is located cannot be determined without further information. Spatial 
frequency is illustrated in Figure 3. 

An extra requirement in the standard DOA problem is a priori knowledge of in- 
coming signal frequencies. Typically, this is handled via a bank of bandpass filters on 
the output of the sensors. Data streams from the sensors at each desired center fre- 
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Figure 3. Spatial Frequency: Two signals with SNR = 2dB, 

0, = 1° and 0 2 = 45° for two snapshots at time = (a) t 0 , ( b ) t v Note 
the variation in 'DC level' as the snapshots sample the nearly broadside 
signal at different phases. 
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quency (f, f , ...) are processed in parallel, resulting in spatial frequencies for each 
temporal frequency. The DOAs are calculated by comparing these spatial frequencies 
with the center frequencies of their respective filter banks. 

Improvement in beamforming may be seen through the use of windows, weighing 
each sensor output by the appropriate amount to narrow the beamwidth or lower 
sidelobes, but, as discussed earlier, at a cost of lowering overall resolution. 

C. SUBSPACE METHODS 
1. Introduction 

The projection-type subspace method utilizes the properties of the 
autocorrelation, covariance, or modified covariance data matrices and their 
eigenvalue/eigenvector decomposition into signal components and noise components in 
estimating the DOA. Generally, subspace methods use an assumed property of the data 
to provide good resolution if the data fits the particular assumptions, even for extremely 
short data sets. If the data (and signals) do not meet the assumptions, the results can 
be quite misleading. The assumptions here call for white noise and a signal whose esti- 
mated autocorrelation matrix converges to the true autocorrelation matrix as the order 
is increased. 

For p complex sinusoids in additive complex white noise the combined 
autocorrelation function of the signal plus noise is given by 



R xx (m) 



p 




( 18 ) 



We define R SJ as the signal autocorrelation matrix of rank p, R„ = a 2 z l of full rank M 
and give the autocorrelation matrix as 



R„ 



p 




(19) 



where I is an M by M identity matrix, R„ is of full rank M due to the noise contribution, 
P, is the power of the /th complex sinusoid, o 2 is the noise variance, and 
e, = [1, e j2nf ‘, ... , eWM-n]? _ Unfortunately, this decomposition is not possible without 
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absolute knowledge of the noise. The p signal vectors e, contain the frequency infor- 
mation and are related to the eigenvectors of R„ by v,. = — ==-e, for i < p. The 

yjM 

eigendecomposition of R„ is 



p 

l 



(/!,• + + 



M 

1 



(20) 



The principal eigenvectors of R„ are a function of both the signal and noise. If no signal 
is present the autocorrelation matrix is simply a diagonal matrix with the eigenvalues 
equal to the variance of the noise [Ref. 1: pp. 422-423]: 



a] 0 
0 a] 0 



0 



0 



. 0 

0 a\ 



(21) 



The data generated from taking M samples of p sinusoids in white noise can be 
used to generate an autocorrelation matrix with the following properties: 

• The theoretical autocorrelation matrix will be composed of a signal autocorrelation 
matrix and a noise autocorrelation matrix. 

• The signal autocorrelation matrix is not full rank if M> p. 

• The p principal eigenvectors of the signal autocorrelation matrix may be used to 
find the sinusoidal frequencies. 

• The p principal eigenvectors of the signal autocorrelation matrix are identical to the 
p principal eigenvectors of the total autocorrelation matrix. 

The matrix will have a minimum eigenvalue ). = a J. [Ref. 1: pp. 422-434] 

Furthermore, the noise subspace eigenvectors are orthogonal to the signal 
eigenvectors, or any linear combination of signal eigenvectors. For the theoretical 
autocorrelation matrix, R u , all eigenvalues not associated with the signals are associated 
with the noise and are identical in magnitude at ). = o] , the minimum eigenvalue of 
R w . Thus, 
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R M V M = 4.I.VM 



( 22 ) 



The zeros of this minimum eigenvector polynomial 



Z 

1=0 



y p+lil+ 0* 7 



(23) 



will have p zeros on the unit circle corresponding to signal frequencies [Ref. 4: pp. 
335-337, Ref. 5: pp. 371-372], 

For the simple case of M - I complex sinusoids, the autocorrelation matrix 
will have a single noise subspace eigenvector \ M with its associated eigenvalue ). = a 2 , , 
the minimum eigenvalue of R jV1 > Thus, the resulting zeros correspond exactly to the 
sinusoidal frequencies. 

2. MUSIC 

Multiple Signal Classification (MUSIC) is a form of a noise subspace fre- 
quency estimator, based on noise subspace eigenvectors with uniform weighting. The 
MUSIC algorithm finds a pseudo spectral estimate from [Ref. 5: p. 373] 



PuusM = 




(24) 



where e, = [1, ... , . The advantage of this method is in its generalized 

nature. There is no longer is a requirement for uniform spatial sampling. Any array 
geometry will provide a solution. A well-designed array will eliminate bearing ambigui- 
ties and provide unique solutions [Ref. 2: pp. 19-28]. 

R.O. Schmidt [Refs. 2, 12. 13] has shown that a group of sensors excited by a 
stationary point source emitter of known frequency will have a magnitude and phase 
relationship (or correspondence) based on the DOA of the plane wave. This corre- 
spondence depends on the array geometry, as well as individual sensors characteristics 
(i.e., directivity, gain, and frequency response), and may not be unique for that one di- 
rection of arrival. 
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By examining the signal in terms of an M dimensional complex field, where each 
sensor provides an orthonormal axis, one can see that a single signal will result in one 
steering vector. This steering vector describes the relationship between the individual 
sensors in terms of phase and magnitude differences, thus l'or any unique signal fre- 
quency and direction of arrival there is one unique steering vector. The magnitude of 
the vector will vary with time, but its direction in M space is a constant determined by 
the amplitude and phase relationship of the sensors for that particular signal as illus- 
trated in Figure 4. 

The theory of superposition applies, thus with two signals present the instanta- 
neous received steering vector is a linear combination of the individual steering vectors. 
The time varying steering vector will move in a plane that is spanned by the individual 
steering vectors. Figure 5 show’s the subspace plane spanned by two steering vectors. 
The same idea can be expanded to a case o f p independent signals causing the steering 
vector to move through a p dimensional subspace (as long as M> p). 

Unfortunately, the steering vector may not determine the actual DOA uniquely. 
In the example of a one-dimensional linear array, a signal gi\es an infinite number of 



Sensor 3 




Figure 4. Steering vector for 3 sensors and 1 signal 
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Sensor 2 



Sensor 1 

Figure 5. Signal Subspace for 2 signals 

DOAs that lie in a cone that is formed by revolving the actual DOA about the axis of 
the array. Thus the array design and its manifold (expected response) will play a large 
part in achieving optimum results. The array manifold describes the predicted response 
of the array to any possible signal or combination of signals. The manif old may be es- 
timated analytically or through physical calibration. 

Analytically describing an array is a complex mathematical procedure for all but 
the simplest arrays (i.e., equally spaced sensors in a linear array or a 3 element 
orthogonal array). It also assumes that absolute knowledge of the array geometry is 
available - an assumption which can lead to errors if the differences are too large. 

Calibration with test sources requires a known, low noise environment while the 
test sources of each desired frequency are placed in a finite number of possible locations. 
The estimated response to actual signals is an interpolation of these responses. Each set 
of calibration parameters requires storage in memory; this results in overall massive 
storage requirements for a comprehensive grid. 

Several parameters such as the number of signals present, the directions of ar- 
rival of those signals, and the signal strengths can be determined from the signal sub- 
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space information. More complex models, however, can be developed that can 
determine signal frequency and polarization as described in References 2 and 13. 

We will now describe the basic steps in the MUSIC algorithm for the DOA 
problem. First, we sample the signals and compute the autocorrelation matrix of the 
data R. Then, we determine the ordered set of eigenvalues and eigenvectors of R. In 
particular, the eigenvectors associated with the minimum eigenvalues must be accurate. 
In the theoretical case the signal eigenvalues are composed of signal strength (P, ) and 
noise variance (a 2 ,) and the nonprincipal eigenvalues will all be identically a 2 v . 

We now determine the number of signals by eigenvalue comparison. A simple 
counting of the eigenvalues greater than a 2 will give the number of signals present in the 
theoretical case. However, the sample autocorrelation estimates does not lead to such 
simple results. Small power level signals may have small eigenvalues (perhaps smaller 
than a 2 ), and the noise eigenvalues will not actually be identical but will group or cluster 
near an approximation of a 2 . Methods of solution include likelihood ratio tests where 
the gaps between the eigenvalues determine threshold placement [Ref. 2: pp. 77-79]. 

Once we find the number of signals present we can determine the signal sub- 
space by the span of the first p eigenvectors 

V/f=[v i v 2 V P ] (25) 

The signal nullspace is its orthogonal complement 

V A -c = [y p+1 v p+2 ... v M ] (26) 

We now find the intersection of the signal subspace with the array manifold. The 
intersection is given in equation 24 for the case of the uniformly spaced linear array. In 
actuality the intersection is estimated with a "best-fit" method used to determine the 
optimum result. In the nonlinear case the array manifold is much more difficult to rep- 
resent. 

Two major areas of difficulty with the MUSIC algorithm are the calculation and 
storage of the array manifold and performing the eigendecomposition of the 
autocorrelation matrix that results from very large arrays. 

3. ESPRIT 

In Reference 14, Paulraj, Roy, and Kailath describe Estimation of Signal 
Parameters via Rotational Invariance Techniques (ESPRIT), a subspace method which 
utilizes two identical subarrays X and Y. A known distance called a displacement vector 
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separates the two parallel subarrays, but no rotation can be present. Each sensor in a 
matching pair (doublet) must be identical, but knowledge of individual sensor and array 
response characteristics is not required. 

The .V elements of both arrays are sampled simultaneously with the signal at 
each sensor being given by 

*/« = 

k=t 



P 

I 



Ski'MQk) + n x,{t) 



= 



z 



,2s£ 

s k {ty >- 



\(6 k ) + n m 



(27) 



where the sampled signal at each sensor in a doublet differs only by a phase term and 
additive noise. With p signals present, s k is the wavefront at the reference sensor in the 
X array, B k is the DOA relative to the displacement vector, a,{d k ) is the response of the 
/th sensor in a subarray relative to the reference sensor in that array for a signal from 
bearing Q k , d is the magnitude of the displacement vector and n the noise term. In vector 
notation the signals at the subarrays are 



x(0 = Hs(/) + n^r) 

\(t) = HcDs(f) + n^/) 

where 

x(0 = [a(0,x 2 (0, ... ,*„ ? (r)] r , 

y(0 = Oi(0.>*2(0 y m i .OH 7 ". 

and ny(r) = L^(t), %('), - » (OH 7 



(28) 



(29) 



The vector of wavefronts at the reference sensor in array X is represented by s(t). All 
displacements and phase differences are based on this sensor. The p by p diagonal ma- 
trix, O, contains the phase delays that occur with each set of doublets 

O = diagfV^ 1 , .... (30) 

where 4> k = ~~ ~ sin 6 k , as shown in equation 13. 
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If is the signal autocorrelation matrix, the subarrav autocorrelation matrix 
is given by 

Rjoc = HR w H r + a\ I (31) 

The cross correlation between subarrays is 

R xy = HR p /> r H r (32) 

where H is the direction matrix whose columns are the direction vectors for the p 
wavefronts 

¥»*) = wej, h 2 (e k ) h„(e k )3 T (33) 

After some manipulation [Ref. 15: pp. 251-253], we obtain the matrix 

r “ “U)" V R ,v 

= HR fr H r -yHR J>() 0 T 'H 7 ' (34) 

= HR W (I - yO r )H r 

In general, there will be p eigenvalues of this matrix. But when y = the /th row 

of (1 — yd? 7 ) becomes zero, leaving p - 1 eigenvalues. Each value of y where this occurs is 
called a generalized eigenvalue (GE). Now, the DOAs can be determined by 



6 k = arcsin(-^j-) (35) 

Due to estimation errors in the calculation of the correlation matrices, some 
error will be induced. Generally, the GE's will be moved off the unit circle and out from 
the origin, but in the case of strong signals, they will still be easily discernible. It should 
be noted that poor array design may result in possible ambiguities (similar to MUSIC). 

Advantages to note over the MUSIC algorithm come with respect to the array 
manifold. With ESPRIT, no manifold is required. The considerable calibration efforts 
and storage requirements are nonexistent. However, the two subarrays must be identical 
in all respects and must be positioned exactly parallel to each other [Ref. 14], 
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4. Other Subspace Methods 

Large variance effects may be seen in the above methods due to the poor reli- 
ability of the estimation of the eigenvectors associated with the minimum eigenvalues 
of the estimated autocorrelation matrix R„. A different method of spectral estimation 
is the use of the principal components where only the eigenvectors associated with the 
largest eigenvalues are used. Some methods have tried to minimize the effect of noise 
by using R„ — a;I but the estimation errors in both R„ and crj have limited their suc- 
cess [Ref. 1: pp. 425]. Other spatial spectral methods include the parametric methods 
such as AR and ML [Ref. 1: pp. 426-427]. 

The structured matrix approximation of Kumaresan and Shaw [Ref. 16] uses K 
snapshots in time of an A' element uniformly spaced linear array with each snapshot in 
time forming a column of a data matrix X, which is then approximated by X M in the least 
squares sense. The bearings information is then calculated using the properties of the 
Vandermonde matrix. 

A combination of frequency domain beamforming and autoregressive modeling 
techniques have been employed in Reference 17 to estimate the direction of arrival in a 
multiple source localization situation using a planar array. 

Halpeny and Childers [Ref. IS] use frequency-wavenumber filters to break the 
multiple wave case down to a succession of single wave problems. 

Reference 19 explains an algorithm that uses non-noise eigenvectors from a 
covariance matrix to obtain high resolution direction of arrival for narrow band sources. 

An adaptive beamforming method similar to a minimum energy approach is 
decribed in Reference 20. The eigenstructure of the correlation matrix is analyzed and 
the computations are modified to optimize resolution but at a cost of array gain. 
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III. THE LANCZOS ALGORITHM 



Lanczos algorithms provide a method to find some of the extreme eigenvalues 
(smallest or largest) and their associated eigenvectors from large matrices with fewer 
operations than required in an entire matrix decomposition. The procedure is a 
tridiagonalization of the original matrix based on an iteration developed by Cornelius 
Lanczos [Ref. 21: pp. 49-206]. Once the matrix is in a tridiagonalized form the 
eigenvalues can be easily computed using a symmetric QR algorithm [Ref. 15: pp. 
278-279] or bisection (with or without Sturm sequencing) [Ref. 15: pp. 305-308], The 
algorithm takes advantage of "minimized iterations" providing quick convergence to the 
final tridiagonal matrix and avoiding accumulation of roundoff error. [Ref. 22] 

The method fell from favor as a tridiagonalization technique with the advent of the 
Givens and Householder transformations later on in the 1950s. Givens transformations 
[Ref. 15: pp. 38-47] use plane rotations (orthogonal matrices) to zero out undesired en- 
tries in the matrix undergoing tridiagonalization. The Householder methods [Ref. 15: 
pp. 43-47] use elementary reflectors to accomplish the same end. Both methods are in- 
herently stable, with the Householder method being slightly superior in both speed and 
accuracy. Both methods outperformed Lanczos as a complete tridiagonalization proce- 
dure in the general case where all eigenvalues are required [Ref. 23: pp. 42-43]. 

The real power of the Lanczos method however lies in obtaining the extreme values 
quickly. The entire matrix need not be tridiagonalized before solutions start to converge. 
Also, if the original matrix is sparse, the Lanczos method maintains that property, sav- 
ing even more computations. Thus for large matrices (order > 100) the Lanczos method 
will converge on extreme eigenvalues in many fewer operations. Recently Tufts and 
Melissinos [Ref. 24: pp. 43-47] have derived a variation of the Lanczos method for high 
resolution spectral estimation and showed that their method outperforms both the sin- 
gular value decomposition and the power method of principal eigenvector and 
eigenvalue computation. Later in this chapter, it will be shown that storage require- 
ments are also minimized. 

This chapter starts with an explanation of the classical Lanczos algorithm as devel- 
oped by Lanczos and modified by Paige [Ref. 25]. Then we will discuss the unorthodox 
Lanczos algorithm of Cullum and Willoughby where no reorthogonalization is per- 
formed [Ref. 23]. Finally, the algorithm used in the direction of arrival problem are 
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discussed in detail. Also, we present some results of the algorithm in the form of spectral 
estimation of multiple tones. 

A. CLASSICAL LANCZOS AND ITS VARIANTS 

The original algorithm u r as designed as a means to directly tridiagonalize the sym- 
metric matrix A. The development of Lanczos algorithm requires the knowledge of 
Krylov sequences and subspaces. For an n by n matrix A and any arbitrary non-null 
n by 1 vector Vj the Krylov sequence of vectors is defined as: 

v i+1 = Av, = A i \ l for k = 1, 2, ... , n (36) 

In the above sequence there will be a vector, say v^, which can be expressed as a 
linear combination of all the preceding vectors. The Krylov matrix of rank m is then 
given by 

V m = [v, v 2 V| ... vj = [v 1? Av„ A 2 Vj, ... , A m-1 v 1 ] (37) 

and the Krylov subspace is the space that spans the vectors [v„ v 2 , ... , vj, 

A^(A,v,) = span Av 1( A 2 v,, ... , A* -i Vj| (38) 

The columns of the n by m Krylov matrix V m are orthogonal. The tridiagonalization 
of the given matrix A is then obtained as 

T=v£av„, (39) 

where T is an m by m tridiagonal matrix. Thus, the given matrix A can be 
tridiagonalized provided we have an efficient way to compute the orthogonal matrix V,, 
or to compute the elements of matrix T by performing the decomposition of equation 
39 directly as proposed by Lanczos [Ref. 22, Ref. 15]. 

The most direct way of performing the tridiagonalization assumes that T = V r AV 
where V = [v 2 v 2 ... v m ] . Note that A is a symmetric nb\ n matrix and V is an 
n by m orthogonal matrix constructed from the Krylov sequence of vectors. The basic 
recursion for a real n by n matrix A starts with a randomly generated unit Lanczos vec- 
tor Vj. Define scalar /?[ = 0 and v 0 = 0. Define Lanczos vectors v, and elements a, and 
/ti for / = 1 , 2 ,..., \I, 



Pi+A'i+i = Av, - « /V/ - /Jp' .i 



(40) 
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o.i = \J Av, (41) 

A+i = v wA (42) 



This is referred to as the basic Lanczos iteration or recursion. The actual Lanczos ma- 
trix T ), j = 1, 2, ... , M, is a symmetric tridiagonal matrix with a, , 1 < / <j, on the di- 
agonal, and /?,*i , 1 < i < {j - 1), above and below the diagonal. 



<*i P 2 0 

V-2 P 3 0 



0 



(43) 



* * • Pj 

0 pj o.j 



Each new Lanczos vector v,^ is obtained from orthogonalizing the vector Av, with 
v, and v,_, . The elements a, and p are the scalar coefficients that make up the Lanczos 
matrix. The basic Lanczos recursion may be condensed into matrix notation by defining 
V„ = ( v i. v 2 * , V J- Then from equations 40, 41, and 42 

AV„ = V„T„ + (44) 

where e m is the coordinate vector whose mth component is 1 while the rest of the com- 
ponents are 0, T m is the m by m Lanczos matrix, the diagonal entries are 

T m (k,k) = a k for 1 < k < m, (45) 

and the entries along the two sub-diagonals on either side of the principal diagonal are 

T m (k.k + 1) = T m {k + 1,/c) = p k+l for 1 < k < m - 1 (46) 

Note that A is never modified (unlike in Givens and Householder transformations), thus 
advantage may be taken of any existing sparsity. Also, the only storage requirements 
are for the three Lanczos vectors (n dimension), the Lanczos matrix T,, and the original 
matrix A. 

We summarize the actual steps: 

• Use the basic recursion of equations 40, 41, and 42 to transform the symmetric 
matrix A into a family of symmetric tridiagonal matrices T, , j= 1,2,..., A/. 
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• For m < M Find the eigenvalues of T m , p (also known as the Ritz values of T m ). 

• Use some or all of these eigenvalues, p, as approximations of the eigenvalues of 

a, . 

• For each eigenvalue p find a unit eigenvector u so that T m u = pu . 

The Ritz vector y is the approximation of the eigenvector of A. It is found from map- 
ping the eigenvector u of the Lanczos matrix. 

3’ 55 v mU (47) 

So the eigendecomposition of T m finally results in the Ritz pair (p, y) which approxi- 
mates the eigenvalues and eigenvectors of A. [Ref. 23: pp. 32-33., Ref. 15: pp. 322-325.] 

Unfortunately, the effects of finite precision arithmetic prevent the theoretical 
Uanczos algorithm from working. The eigenvalues of A and the eigenvalues of T m no 
longer converge. Roundoff errors are partially to blame, but the dominant effect is from 
the loss of orthogonality of the Uanczos vectors v, . From equation 40 it can be seen 
that a small will have great effect on v,^. Paige showed that the algorithm runs 
within allowable error constraints until it starts to converge on the first eigenvalue. At 
this point becomes small and the Uanczos vectors lose orthogonality. The loss of 
orthogonality is not random, however. It always occurs in the direction of the Ritz 
vector y. 

This trouble can be dealt with through reorthogonalization. But questions that we 
must answer include: 

• How much reorthogonalization is required? 

• Where should it be performed? 

• Reorthogonalize with respect to what? 

Complete reorthogonalization is the first choice, inserting a Householder matrix 
computation into the Lanczos algorithm. This enforces orthogonality among all the 
Lanczos vectors and is effective at keeping the system stable. But the extra computa- 
tions it requires negate any advantage of the Lanczos algorithm, making the number of 
computations on the same order as a complete decomposition. Numerous vectors have 
to be stored and compared requiring many swaps into and out of storage. [Ref. 15: pp. 
334-335] 
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Selective reorthogonalization is clearly more efficient. The extra computations are 
performed only if orthogonality checks indicate loss is imminent. Paige shows that the 
loss of orthogonality occurs only when the algorithm converges on a Ritz pair. Thus, 
instead of reorthogonalizing against every other Lanczos vector, using the less numerous 
Ritz vectors will be sufficient. Storage is therefore minimized. Only when all eigenvalues 
of the matrix are required does this method become computationally more expensive 
than other techniques. [Ref. 26] 

The last option is no reorthogonalization. The explanation above would seem to in- 
dicate that maintaining orthogonality is required. However by analyzing the causes and 
effects of the original loss of orthogonalization one can insert corrections into the sol- 
ution to give valid eigenvalues and eigenvectors. This is the procedure that will be ex- 
amined in the next section. 

One other property will be mentioned. Here the single vector Lanczos recursion 
described above will not find duplicate eigenvalues of the matrix A. Parlett's proof uses 
the power method to expand v to compute the columns of the Krylov matrix K”(v, A). 
If J{\) is the smallest invariant subspace of [R -1 which contains v, then the projection of 
A onto S(x) has simple eigenvalues. Also, the Krylov subspaces fill up J(x) and for 
some l < n we have 

span{x) cz A 2 ( A. v) cr A' 3 (A, v) S ... m K l ( A, v) = A' /+1 (A, v) = J(x) (4S) 

The result is that some eigenvectors of A may be missed, and any repeated eigenvalues 
will not be detected. [Ref. 27: pp. 235-239] 

The multiple eigenvalue problem can be treated by using a Block Lanczos algorithm. 
Instead of obtaining a tridiagonal matrix, the result is a banded block matrix, where the 
diagonal is M„, an / by / matrix, and the above and below r the principal diagonal are 
upper triangular matrices B„ . Each block must be dimensioned / > p , where p is the 
estimated multiplicity of the desired eigenvalue. 



Mi B 2 r 0 

B 2 M 2 B 3 r 0 

0 B 3 . . . 



(49) 



0 

0 
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This is analogous to the general case [Ref. 15: pp. 337-339]. The block Lanczos algo- 
rithm is briefly reviewed at the end of this chapter. 

The above discussion assumes that the given matrix A is real and symmetric. Besides 
the algorithms summarized in this section for a real and symmetric matrix case, there 
are other general Lanczos algorithms proposed in the literature [Ref. 23] that are suitable 
for Hermitian matrices, non-symmetric matrices, and for rectangular matrices. 

B. IMPLEMENTATION 

The Lanczos phenomenon states that a few, not all, of the eigenvalues of a very 
large matrix A can be computed using the Lanczos recursion with no 
reorthogonalization. But to find most of the eigenvalues, the Lanczos matrix, T m , will 
grow in size approaching that of A, causing the loss of orthogonality of the Lanczos 
vectors. The loss of orthogonality results in many spurious eigenvalues, as well as extra 
copies of good eigenvalues. In any case, a test is required to confirm either: 

• a "found" eigenvalue is good, or 

• the eigenvalue that appears is spurious. 

Golub and Van Loan, Parlett, and Paige [Refs. 15, 27, and 28] describe procedures 
that look at the eigenvalues for each T m as m is stepped up in size. All the eigenvalues 
of T m are computed at each step. The good eigenvalues will repeat at each larger T m , 
while the spurious eigenvalues jump around. If an eigenvalue does not match at con- 
secutive T m 's it may be considered spurious and thrown out. If a good eigenvalue is 
mistakenly deleted (due to numerical roundofl), it can be counted on to reappear in the 
next step. 

Cullum and Willoughby [Ref. 23] take a different tack by developing a test to find 
and eliminate bad eigenvalues and retain the rest. The advantage here is in utilizing the 
machine's tolerance to drop bad eigenvalues, while not discarding good eigenvalues that 
have yet to converge. As a result a larger spectrum is available sooner, even though it 
may only be a rough estimate of where the eigenvalues will finally converge. 

In practice, parts of the Lanczos recursion (equations 41 and 42) are replaced by 

a , = - Avm) (50) 

and 

= ||AVi — (51) 
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Computation of the element a, is a modified Gram-Schmidt orthogonalization proce- 
dure. The new /? (+ , is equivalent to the /?, +1 of equation 42 but now it directly controls 
the size of the Lanczos vector. 

In what follows we describe two Lanczos algorithms, namely the single vector 
Lanczos algorithm which is modified and analyzed by Paige [Ref. 28] and the block 
Lanczos algorithm described by Cullum and Willoughby [Ref. 23]. Both algorithms have 
been considered for the estimation of the directions-of-arrival of multiple targets in noisy 
environments in this thesis. 

1. Single Vector Lanczos Algorithm 

The first procedure to be described is the Paige's single vector Lanczos algo- 
rithm [Ref. 28] for real symmetric matrices. The single vector Lanczos procedure is one 
of the most straightforward implementations of the theory. This procedure will find 
some or many of the eigenvalues and eigenvectors of a real symmetric matrix A such that 
Ax = Ax. It will not detect repeated eigenvalues. However, it may be noted that for 
many problems of interest in practice we do not have strictly multiple eigenvalues. For 
example, in the direction-of-arrival problems the smallest eigenvalues of the 
autocorrelation matrix corresponding to the noise associated with the target signals are 
spread over a small range rather than coinciding on the same value [Ref. 2]. 

No reorthogonalization is performed as part of the single vector Lanczos al- 
gorithm. As mentioned earlier, the Lanczos vectors begin to lose their orthogonality 
when we seek to estimate all or most of the eigenvalues of the real symmetric matrix A. 
For the application under consideration, however, we are generally interested in only a 
few of the minimum eigenvalues and the corresponding eigenvectors. It is mainly for this 
reason that we have not attempted the complete or selective reorthogonalization of the 
Lanczos vectors in this study. 

Now we shall outline the basic steps involved in the single vector Lanczos al- 
gorithm. This is based on the recursion described by equations 40, 50 and 51. Based on 
these equations Paige [Ref. 28] presented four different single vector algorithms. We 
have adapted one of in this study. The complete eigenvalue/eigenvector problem actually 
consists of three parts: (a) obtaining the tridiagonal matrix T m from the given symmetric 
and real matrix A using Paige's recursion, (b) determining the smallest eigenvalues of 
T m using the bisection method and Sturm sequencing, and (c) estimating the corre- 
sponding eigenvectors by computing the Ritz vectors. The details are presented in the 
following. 
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Step 1: As shown in equation 43 the symmetric tridiagonal matrix T, has entries 
<Xj and Pj along its principal, and the adjacent sub and super diagonals, respectively. The 
following recursive expressions are then used to compute the entries of the tridiagonal 
matrix, and also the Lanczos vectors v, [Ref. 28]: 

Initial conditions: v, is an arbitrary/? by 1 vector such that ||v|| 2 = 1 



where w and u r are some intermediate vectors. The vector v, is obtained by filling its 
entries with a random number sequence and then normalizing it with respect to its 
Euclidean norm. Now T m is obtained by simply filling its entries as in equation 43. Note 
that m <£ n in the above. One quick test to ensure that we have obtained a fairly accurate 
estimate of T m is to compute the product v/vj . The result should be equal to b H , where 
S, ; is the Kronecker delta function. 

Step 2: The eigenvalues of the tridiagonal matrix T m , denoted by fx p can be 
computed using the bisection method and Sturm sequencing. Actually one could obtain 
both eigenvalues and eigenvectors of T r by employing such methods as the QR algo- 
rithm. However, when only a few eigenvalues are required, the bisection method seems 
appropriate. For the given m by m matrix T m we define the characterstic polynomial 



u, = A Vl 
<*o = h = 0 



for j = 1, 2, ... , m 




(52) 



"J = u y - 

Pj = lh>,n 2 

\ J+l = ysjlPj 

u /+> = Av m - fy'j 



PM = det(T w - M I) 



(53) 



which can be recursively computed as follows 



Po(p) = 1 
P\{p) = «i ~P 

PM) = i a j — P)pj-\(P) ~ P]-\Pj-2(P-) 



(54) 



for j = 2, 3, ... , m 
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The roots of the polynomial p m (p) are the required eigenvalues. For our appli- 
cation, we are only interested in a small range of eigenvalues at the lowest end of the 
eigenvalue spectrum. We make use of the Sturm sequencing property that the 
eigenvalues of T,._[ strictly separate those of T ; [Ref. 15: pp. 305-307] and implement the 
following iteration: 

PL - a + b 

b = P if pjfl)p m {b)< 0 (55) 

a = P if Pm( a )Pm(b) > o 

and we repeat the above as long as | b — a | > e( | b \ + | a \ ), where c is the machine 
unit round-off error and \_a,b ] is the range of our required eigenvalues. Determining the 
range of interest in our application may require some a priori knowledge about the signal 
to noise ratios (SNR) and it may take a couple of iterations to do this. Some alternatives 
to the iteration given in equation 55 are to use a straightforward polynomial root finder 
and then pick the roots of interest, or to employ the well known L-D-U factorization, 
both of which may not be computationally efficient. 

Step 3: There are two ways to obtain the eigenvectors of A knowing its 
eigenvalues, ). r Note that p ; are the estimates of J| In the first method, we compute 
the eigenvectors of T m , denoted by t„ and then obtain the eigenvectors of A given by 

Xj = L m t| (56) 

where L m = fv, v 2 ... v m ] is the Lanczos matrix which ideally is the same as the 
Krylov matrix of equation 37. Note that {p,, t y .) e T m . 

The second method involves computing the Ritz vectors either from the 
Rayleigh quotient interation or by the orthogonal iteration. Here we assume that we 
have good estimates of X t from the previous step, and proceed to obtain the eigenvector 
Xj by minimizing the cost function 



J = ||(A - > v I)x y || 2 (57) 

It can be shown that a simple minimization of J with respect to x, yields the Rayleigh 
quotient of x, 
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r(xj) = Ij = 



(58) 



x‘ Axj 



Therefore, given )., and using equation 58 we can formulate the Rayleigh quo- 
tient iteration as follows [Ref. 15, 27] 

Initial condition: Xq is an arbitrary vector such that ||x 0 || 2 = 1 
for k = 0, 1,2,... 



r(x ; ) = 



x k Ax k 



(59) 



solve (A-r(x fc )I)y* +1 = x k for y k+1 



*k+{ = y*+i/lly*+ill2 

where y A _, is some intermediate vector. We stop the iteration when ||y^,|| 2 converges to 
a constant or when r(xj ^ /*, one of the known eigenvalues. At each iteration step we 
need to solve an n by n system of equations in this method. One advantage with this 
method, however, is that it converges very quickly. Besides the above iteration, some 
other methods are outlined in References 15 , 23 and 27. We remark that if only a few 
eigenvalues and eigenvectors (say, five) are required, it may be more direct to use 
equation 56. 

We now present an example of the ability of the single vector Lanczos algorithm 
to estimate the direction of arrival, or to find spectral lines in noise, and the advantages 
in extracting more than one eigenvalue and eigenvector in this process. We consider a 
signal with three sinusoids present in noise 



*(») = 



3 




A ( cos(27r/J«) + n(n) 



(60) 



where A, are the amplitudes of sinusoids, f are the normalized spatial frequencies 
(0 <f < 0.5 corresponding to .0 < d < ) , and n{n) is the zero mean white noise with 

a variance of a]. 
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We have computed a 25 by 25 autocorrelation matrix of *(«), R„ , by using 100 
data samples. We have used the covariance method for this purpose, hence R„ is real 
and symmetric. The eigenvectors x y of R„ corresponding to the lowest eigenvalues are 
computed by employing the single vector Lanczos algorithm. The power spectral density 
estimates are computed as follows: 



where x Jt are the elements of the jth eigenvector, x y . 

Figure 6 and Figure 7 show the power spectral density (PSD) estimates of 
equation 60 with an SNR of 10 dB for j = 1 and 2, respectively. 
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(61) 
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Figure 6. PSD of first eigenvector 



Note that the index j indicates the increasing magnitude of the eigenvalues. Thus, 
(/l„x,) are the lowest eigenvalue and the corresponding eigenvector. In both figures we 
have the peaks at the correct locations (9°, 27°, 63°). However, they both have spurious 



smaller peaks at different locations, We can observe the same trend for the first five 
eigenvectors as shown in Figure 8, where the spectral estimates are over laid on each 
other. 

Based on the above results one feels that we could employ some kind of aver- 
aging to get rid of the spurious peaks and improve the estimation performance. We have 
implemented two such methods: the eigenvector averaging and the spectral averaging. 
Figure 9 shows the result of the algebraic averaging of the first 5 eigenvectors, and 
Figure 10 shows the result of the algebraic averaging of the spectral estimates of the 
same eigenvectors. As seen from Figures 9 and 10, eigenvector averaging yields better 
results than spectral averaging. 




Figure 7. PSD of second eigenvector 

Improved results, however, were obtained by using what is called spectral 
multiplication which is obtained by taking the product of the individual spectra, given 
by 

j 

s«W = []sj£w («) 

J= 1 
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where 7 is a predetermined number (J <m< n). Figures 1 1 and 12 show’ the results of 
spectral multiplication for 7 = 2 and 7=5, respectively. As can be seen in these fig- 
ures, using more spectra in equation 62 greatly improves the performance. Also, even 
for 7=2, spectral multiplication outperforms the eigenvector averaging method. 
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Figure 8. Overlayed PSDs of first 5 eigenvector 

In the remainder of the thesis, we have used spectral multiplication in prefer- 
ence to the eigenvector or spectral averaging. Figure 13 shows the multiplication of 5 
spectra for the case when the SNR = 0 dB. We notice a spurious peak around 
6 — 74°. More spurious peaks are observed when the SNR is decreased to -5 dB (see 
Figure 14) and Figure 15 shows the spectrum for the SNR but we have used the 
eigenvectors 6-10 in this case. Improved performance is obtained as shown in 
Figure 16 (J = 10) and Figure 17 (J = 15) by using more and more eigenvectors in the 
spectral multiplication. 

In all the above cases we always observed the signal spectral peaks at the right 
places. The spurious peaks, however, did not appear at the same location as we used a 
different eigenvector to compute the spectrum, S^’(/). 
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Figure 9. Eigenvector averaging 




Figure 10. Spectral averaging 
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Figure 11. Spectral product for 2 eigenvectors 




Figure 12. Spectral product for 5 eigenvectors 
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Figure 13. Spectral product for 5 eigenvectors, 0 dB 




Figure 14. Spectral product for 5 eigenvectors, -5 dB: Using second through sixth 

eigenvectors 
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Figure 15. Spectral product for 5 eigemectors, -5 dB: Using sixth through tenth 

eigenvectors 




Bearings 

Figure 16. Spectral product for 10 eigenvectors, -5 dB 
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2. Other Methods 

The single vector Lanczos algorithm will not determine that repeating 
eigenvalues exist, thus it cannot find the corresponding eigenvectors. The subspacc that 
results has an incomplete basis as it is described only by the eigenvectors that are com- 
puted. The solution is to use a block method that is analogous to the single vector 
Lanczos algorithm. As we mentioned earlier, the block form of the Lanczos algorithm 
does fmd eigenvectors with multiplicity p as long as the blocks are dimensioned / ^ p. 
We attempted to incorporate the Cullum and Willoughby hybrid block Lanczos algo- 
rithm (Ref. 29 Chapter 8) into our direction of arrival model. W'e postulated that it 
would be desirable to compute a few of the extreme smallest repeating eigenvalues and 
their respective eigenvectors. However we were never able to get the program to reliably 
compute good eigenvalues and eigenvectors for the autocorrelation matrix. This has not 
posed a problem for our model as the covariance matrix does not appear to have re- 
peating eigenvalues, but a larger order matrix may indeed include duplicating noise 
eigenvalues and require an algorithm that will accurately operate with that perturbation. 
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The algorithm we attempted to use is actually a hybrid approach to finding the 
eigenvalues and eigenvectors of a real symmetric matrix A. For insight into the problem 
look at the block analogy of equations 40, 41, and 42. Define matrices 
B! = 0 and V 0 = 0. The n by q matrix Vj has columns that are orthonormal random 
vectors. The value of q must be greater than or equal to the number of eigenvalues to 
be found. 



for i = 2, 3, ... , s 

t (63) 

V W B W - P, = AV,-V,M,-V ( _,Br 

M, = V, r (AV,-V i _ 1 B, r ) (64) 

V /+ ,B i+1 S P, (65) 

The matrix B,_, is a modified Gram-Schmidt orthonormalization of the columns of P,. 
Also note that the matrix M, correspond to the «'s of the single vector Lanczos. 

The block analogy to the Krylov subspace approach can be performed with 

K% A.Vj) = sp««{v lf AV lt A 2 V„ ... , A J_1 V} (66) 

The blocks V,, for j = 1, 2, ... , s form the orthonormal basis of the Krylov subspace. 

It can be shown that for a symmetric nbyn matrix A and an orthonormal 
n by q starting matrix V,, that the block recursion equations 63, 64, and 65 will generate 
blocks V 2 , V 3 , ... , \ s where qs < n. It is these blocks that form an orthonormal basis 
of the subspace A' : (A, V,) . In much the same way as the single vector Lanczos algorithm 
generates the tridiagonal Lanczos matrix, the block variant generates blocks, but these 
are now nontridiagonal. At the end of each iteration the Lanczos matrix is of the form 



40 




m[am 



b: 



a 2 @3 



P 3 a 3 P$ 



T s = M 7 AM 



P 4 a 4 



( 67 ) 



• Pt 



The submatrix MfAM consists of the reorthogonalized terms and M, is the portion of 
the first block that is not generating descendants. Ritz vectors are computed on every 
iteration and are used as the starting blocks for the next iteration. Each block is required 
to be reorthogonalized with respect to the all the vectors in first block which is not being 
allowed to generate descendants. It is apparent that the block procedure requires a great 
amount of storage and is very computationally intensive. 
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IV. RESULTS 



Using the Lanczos algorithm it is possible to find some of the eigenvalues and 
eigenvectors of a matrix without going through an entire matrix decomposition. The 
smallest eigenvalue of the autocorrelation matrix and its corresponding eigenvector will 
have the required spectral information to determine a source's bearing (direction of ar- 
rival) from an array. Multiplying several of the resultant eigenvectors' power spectral 
densities will tend to reinforce the true spectral peak and zero out spurious peaks that 
do not occur with every eigenvector. 

The problem with finding the split between the noise and signal eigenvalues disap- 
pears as only a few of the smallest eigenvalues of a large matrix (in relation to the 
number of sources) are used. 

A. APPROACH 

The received signal is modeled by sinusoids at normalized spatial frequencies pro- 
portional to their bearings from endfire (.0 = 0° , .5 = 90° ). The sum of these 
sinusoids is sampled at a rate based on the interelement spacing of 2. Thus a source 
at endfire is sampled at the Nyquist rate and the sample rate increases as the bearing 
shifts toward the array broadside. The simulation uses 



where ss(n) is the instaneous excitation for the sensor at location n, A is the amplitude 
of each of the T signals, f is the normalized spatial frequency of the /th source (de- 
pendent upon bearing), and n[n) is white gaussian noise. The relationship between A 
and the noise variance o\ is determined by the desired signal to noise ratio (SNR), w T here 



The experiment consists of simulating a linear array with equally spaced sensors re- 
ceiving signals of known temporal frequency from various bearings. One possible 



T 




( 68 ) 




(69) 
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physical implementation would place a bank of bandpass filters on each sensor with the 
outputs from each similar filter tied into a correlator. Advantages of this method include 
the processing gain found by prefiltering the noise and simple parallel implementation 
with separate channels for each frequency band. The lowering of the noise bandwidth 
will raise the SXR at the correlator. As more filters are used (smaller bandwidth) the 
noise power decreases and the SXR is increased. The algorithm creates an 
autocorrelation matrix with the output of the correlator. The Lanczos tridiagonalization 
and eigendecomposition provides the eigenvectors that are estimates of the spatial PSD. 
The PSD corresponds to the sources directions of arrival. A possible implementation is 
shown in Figure 18. 

B. EXPERIMENT SET UP 

The first three cases show the effect of different signal strengths on the ability to 
accurately determine the number of targets and the bearing resolution for various di- 
rections and target spacing. In each of these cases, the number of sensors is 100, a ma- 
trix size of 25 is used and 15 iterations (the number of eignevalues eigenvectors found) 
are performed. Case 4 uses three 5 dB sources at 18°, 36° and 41.4° to illustrate the ef- 
fects on changing the number of sensors (samples), the size of the autocorrelation ma- 
trix, and the number of eigenvectors used. The noise is randomly generated white 
gaussian noise with a standard deviation selected to provide the desired SXR. Each 
figure shows, (a) the PSD of selected eigenvectors overlayed and plotted versus bearing 
and (b) the product of selected PSDs of those eigenvectors. 

Case 1 is with all sources at a signal strength of 5 dB. Figure 19 shows results from 
the second through sixth eigenvectors of a single source at 18°. Note that some of the 
eigenvectors have individual peaks as high as the true signal peak, but only at the true 
bearing do all have a common peak. Figure 20 illustrates the other end of the spectrum, 
at 81°. Once again the second through sixth eigenvectors are overlayed to show that the 
correct bearing is consistently displayed, but in this case one eigenvector has an indi- 
vidual peak higher than the true signal peak. The product of these PSDs provides suf- 
ficient resolution. Figure 21 has two closely spaced sources at 36° and 38°. Resolution 
is achieved but the PSD product of the second through sixth eigenvectors shows a spu- 
rious peak near 75°. Figure 22 is from three sources at 0°, 36° and 88,2°. The individual 
eigenvector PSDs clearly show the excellent performance at broadside. 

Case 2 lowers the signal strength of all sources to 0 dB. Figure 23 shows results 
from the second through sixth eigenvectors of a single source at 18°. Many more peaks 



43 




BEARINGS 



Figure 18. A physical implementation 

are visible in the PSD product, making the decision of how many targets more difficult. 
Figure 24 shows that at the the other end of the spectrum (at 81°) the situation is 
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Figure 19. Case 1 5 dB, 1 target at 18 (a) overlay of PSDs from second 

through sixth eigenvectors, (b) product of the above I’SDs 
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Figure 20. Case 1 5 dB, 1 target at 81 °s (a) overlay of PSDs from second 
through sixth eigenvectors, (b) product of the above I’SDs 
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Figure 21. Case 1 5 dB, 2 targets at 36 0 and 38 °: (a) overlay of PSDs from third 
through seventh eigenvectors, (b) product of the above PSDs 
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Tigure 22. Case 1 5 tlli, 3 targets at 0 °, 36 ° and 88.2 (a) ovcilay of PSDs 

from second through sixth eigenvectors, (b) product of the above PSDs 
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slightly worse (due to a lower sampling rate). The overlays of the second through sixth 
eigenvectors show that the correct bearing is consistently displayed, but in this case 
enough spurious peaks reinforce one another, resulting in the PSD product that has not 
zeroed out the bad peaks. Figure 25 illustrates that the proper choice of eigenvectors 
will resolve this problem. Here the PSDs of the first through fifth eigenvectors are used, 
giving a product that is easier to determine correctly. Figure 26 shows the 0 dB case for 
two close spaced sources at 36° and 38°. Resolution is achieved but the PSD product 
of the first through fifth eigenvectors shows several spurious peaks, including the same 
one as in Case 1 near 75°. Figure 27 is the three source example at 0 dB. The individual 
eigenvector PSDs are repeating at the proper bearings but the performance at broadside 
is resulting in the product at the other bearings actually being driven down. 

A signal strength of -5 dB is used for Case 3. Figure 28 shows results from the first 
10 eigenvectors of a single source at 18°. Using more good eigenvectors increases the 
likelihood that all spurious peaks will be diminished. Figure 29 illustrates results at the 
other end of the spectrum (at 81°). Resolution is looked at in Figure 30. At -5 dB the 
algorithm cannot separate the two close spaced sources at 36° and 38°. A number of 
spurious peaks are higher than the bump at 36° making it impossible to accurately de- 
termine the number of sources as well as both locations. Resolution is tried again in 
Figure 31 with 2 sources at 36° and 40° . Using five eigenvector PSD products produces 
good results. Figure 32 shows 3 sources at -5 dB. Good performance is seen both in 
the overlays and in the PSD product. 

Case 4 starts with 100 sensors and a 25 by 25 covariance matrix shown in 
Figure 33. As the number of sensors is decreased and the number of eigenvectors used 
is held constant, more spurious peaks start to occur (Figure 34 and Figure 35). 
Figure 36 shows the spectral improvement as more eigenvectors are used. As the 
number of sensors is decreased to 40 (Figure 37) and seven eigenvectors are used, the 
results are still acceptable. Using only 30 sensors we can no longer see resolve between 
the two closely spaced sources. With a sufficient number of eigenvectors no spurious 
peaks are present but the true number of targets is nondeterminable (Figure 38 and 
Figure 39). 
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Figure 23. Case 2 0 dB, 1 target at 18 °: (a) overlay of PSDs from second 
through sixth eigenvectors, (b) product of the above PSDs 
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Figure 24. Case 2 0 dB, 1 target at 81 °: (a) overlay of PSDs from first through 
fifth eigenvectors, (b) product of the second through sixth eigenvector 
PSDs 
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Figure 25. Case 2 0 dB, 1 target at 81 °: product of the first through fifth 

eigenvector I’SDs 
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Figure 26. Case 2 U df3, 2 targets at 36 ° and 38 °: (a) overlay of FSDs from first 
through third eigenvectors, (b) product of the first through fifth 
eigenvector PSDs 
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Figure 27. Case 2 0 dB, 3 targets at 0 °, 36 0 and 88.2 °: (a) overlay of PSDs 
from second through sixth eigenvectors, (b) product of the above PSDs 
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Figure 28. Case 3 -5 dI3, 1 target at 18 (a) overlay of PSDs from first through 

sixth eigenvectors, (b) product of the first 10 eigenvector PSDs 
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Figure 29. Case 3 -5 dB, 1 target at 81 (a) overlay of PSDs from first through 

sixth eigenvectors, (b) product of the first 10 eigenvector PSDs 
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Figure 30. Case 3 -5 dB, 2 targets at 36 0 and 38 (a) overlay of PSDs from first 

through sixth eigenvectors, (b) product of the above PSDs 
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Figure 31. Case 3 -5 dB, 2 targets at 36 0 and 40 (a) overlay of PSDs from First 

through sixth eigenvectors, (b) product of the first 5 eigenvector PSDs 



58 




Figure 32. Case 3 -5 dB, 3 targets at 0 °, 36 0 and 88.2 (a) overlay of PSDs 

from first through sixth eigenvectors, (b) product of the first 10 PSDs 
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Figure 33. Case 4 5 dB, 3 targets at 18 36 0 and 41.4 °: 100 sensors, (a) PSDs 

of second through sixth eigenvectors, (b) product of the above PSDs 
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Figure 3-4. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 °: 75 sensors, (a) PSDs 
of second through sixth eigenvectors, (b) product of the above PSDs 
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Figure 35. Case 4 5 dB, 3 targets at 18 °, 36 0 and 41.4 °: 50 sensors, (a) PSDs 
of second through sixth eigenvectors, (b) product of the above PSDs 
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f igure 36. Case 4 5 (IB. 3 targets at 18 ", 36 0 and 41.4 °: Product of the first 

through eighth eigenvector PSDs 
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Figure 37. Case 4 5 dB, 3 targets at 18 °, 36 ° and 41.4 °: 40 sensors, 20 by 20 
matrix, (a)PSDs of first through seventh eigenvectors, (b) product of 
the above PSDs 
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Figure 38. Case 4 5 dB, 3 targets at 18 °, 36 0 and 41.4 °: 30 sensors, 20 by 20 

matrix, (a)PSDs of first through fifth eigenvectors, (b)Product of the 
first through fourth eigenvector PSDs 
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Figure 39. Case 4 5 dB, 3 targets at 18 °, 36 0 and 41.4 °: 30 sensors, 20 by 20 
matrix, (a)Product of the first through fifth eigenvector PSDs, 
(b)Product of the first through twelfth eigenvector PSDs 
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v. CONCLUSIONS AND RECOMMENDATIONS 



The results plotted in Chapter IV indicate that the eigenvectors found using the 
Lanczos algorithm are sufficiently accurate to determine the spectrum. Although no 
direct comparisons with other eigendecomposition methods are performed, the theory 
indicates that many fewer operations are required. We handle the other difficulty of 
conventional subspace methods by using only a few of the eigenvectors associated with 
the minimum eigenvalues of the autocorrelation matrix. No estimation of the noise 
subspace dimension is required or performed. 

This theory may be applied to any system requiring rapid decomposition of the 
correlation matrix. Examples include phased array radar and passive acoustic arrays 
[Refs. 30, 31]. Reference 32 details an experimental system using the MUSIC algorithm 
for multiple source direction finding. 

The following areas are recommended for future study. 

• Use of the products of multiple spectra apparently resulted in good detection at low 
SNR. More research in this area to determine a physical interpretation of this 
method is required. 

• Analysis and comparison of the results in terms of computational speed and accu- 
racy with other eigendecomposition methods should be performed to find the true 
results. 

• The Lanczos algorithm developed uses no reorthogonalization nor will it find re- 
peating eigenvalues. Other forms of the Lanczos algorithm are available. Com- 
parisons between these different methods to determine accuracy and speed may 
lead to more optimum results. 

• A more detailed model should be developed that will simulate an array with a bank 
of bandpass filters to better forecast results of a physical implementation (as seen 
in Figure IS of Chapter IV). 

• A method which implements the algorithm in parallel fashion may be tried. Using 
one long linear array, several overlapping subarrays may be used to simultaneously 
create several autocorrelation matrices. The algorithm may be applied to these 
matrices in parallel. It is predicted that the greater number of available 
eigenvectors will more properly describe the noise subspace and therefore more 
accurately estimate the spectrum. 
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